\chapter{谱分析}
\section{引言}
	一个余弦曲线可以书写如下，
	\begin{equation}\label{cos}
	 R\cos (2\pi ft+\Phi) 
	\end{equation}

	
	其中$ R $为振幅，$ f $是频率，$ 1/f $是周期，$ \Phi $是相位。
	
	问题是对于一个组合曲线，如
	\[ Y_t=2\cos \left(2\pi t \frac{4}{96}\right) +3\cos \left[2\pi\left(\frac{14}{96}+0.3\right)\right]\]
	又如何估计其周期性呢？我们知道可用三角恒等式将\eqref{cos}式重新参数化为，
	\[  R\cos (2\pi ft+\Phi) =A\cos(2\pi f t)+B\sin (2\pi ft) \]
	
	那么对于已知的频率$ f $用OLS就可以估计出$ A\text{和}B $。实际上，对于由任意振幅、频率和相位的余弦曲线的一般线性组合可以写成，
	\begin{equation}\label{1-2}
	A_0+\sum_{j=1}^{m}[A_j\cos(2\pi f_j t)+B_j\sin (2\pi f_jt)]
	\end{equation}
	
	OLS就可以估计出$ A $和$ B $。特别地，若观测个数$ n $是奇数，则形如$ \frac{1}{n},\frac{2}{n},\cdots,\frac{k}{n}(=\frac{1}{2}-\frac{1}{2n}) $的频率称为\textbf{傅立叶频率}，这些频率上的余弦和正弦是正交的，因此$ A\text{和}B $的OLS估计为，
	\begin{align*}
	\hat A_0&=\bar Y\\
	\hat A_j&=\frac{2}{n}\sum_{t=1}^nY_t\cos (2\pi tj/n)\\
	\hat B_j&=\frac{2}{n}\sum_{t=1}^nY_t\sin (2\pi tj/n)
	\end{align*}
	
	若样本量为偶数，上述估计对$ j=1,2,\cdots,k-1 $时仍成立，只是在$ j=k $时，有，
	\begin{align*}
		\hat A_k&=\frac{1}{n}\sum_{t=1}^n(-1)^tY_t\\
		\hat B_k&=0
	\end{align*}
	
	要注意，只要我选取足够多的$ m $（即自变量），那么对于任意长度的观测值，我总能完美拟合\eqref{1-2}式。
	\section{周期图}
	周期图如下定义，当样本量为奇数，
	\[ I(f)=I(\frac{j}{n})=\frac{n}{2}(\hat A_j^2+\hat B_j^2) \]
	
	样本量为偶数，注意修改在极端频率$ f=k/n=1/2 $上，
	\[ I(\frac{1}{2})=n\hat A_k^2 \]
	
一般绘制横轴为频率，纵轴为周期图的图形。周期图的好处在于，因周期图与相应振幅的平方和成正比，周期图的高峰就显示了不同频率行余弦-正弦的相对强度。如下图所示，
\begin{figure}[H]
\includegraphics[scale=0.8]{13-1.png}
\end{figure}

注意观察，周期图的横轴区间仅仅限于$ [0,0.5] $，这是因为$ [0,0.5] $上的每一个频率与形如$ f+k\cdot 0.5 $的频率互为假频，即在离散的时间点上，看不出它们图形的区别。

另外，尽管周期图是在离散的傅立叶频率上计算的，实际上可以将$ f $在$ [0,0.5] $的区间上连续化。
\section{谱表示和谱分布}
\subsection{谱密度}
考虑一个一般的时间序列，总可以写成，
\begin{equation}\label{3-1}
Y_t=\sum_{j=1}^{m}[A_j\cos(2\pi f_j t)+B_j\sin (2\pi f_jt)]
\end{equation}

其中，$0< f_1<f_2\cdots <f_m<1/2 $固定，$ A_j,B_j $是独立正态随机变量，均值为0，$ Var(A_j)=Var(B_j)=\sigma_j^2 $。可知，$ Y_t $均值为0，且平稳，自协方差为，
\[ \gamma_k=\sum_{j=1}^m\sigma_j^2\cos(2\pi kf_j) \]

\textbf{谱分析，习惯上要先移去样本均值}。定义，在$ [-\frac{1}{2},\frac{1}{2}] $上的频率，\textbf{样本谱密度}为$ \hat S(f)=\frac{1}{2}I(f) $，且$ \hat S(1/2)=I(1/2) $，经过复杂的代数运算，可以有，
\[ \hat S(f)=\hat \gamma_0+2\sum_{k=1}^{n-1}\hat \gamma_k\cos (2\pi fk) \]

可见，样本谱密度是样本协方差的傅立叶变换。也可以看到样本谱密度与样本协方差是一一映射的，它们包含的信息相同。而且，样本谱密度曲线下方的面积即是时间序列的样本方差，$ \hat \gamma_0=\int_{-1/2}^{1/2}\hat S(f)df=\frac{1}{n}\sum_{t=1}^n(Y_t-\bar Y)^2 $。实际上，理论谱密度可以写为，
\begin{equation}\label{3-3}
S(f)=\hat \gamma_0+2\sum_{k=1}^{n-1}\gamma_k\cos (2\pi fk) 
\end{equation}
\subsection{线性滤波器}
对时间序列$ \{X_t\} $滤波产生$ \{Y_t\} $，即，
\[ Y_t=\sum_{j=-\infty}^{\infty}c_jX_{t-j} \]

实际上，AR，MA，ARMA过程都是对白噪声的一个线性滤波。上式通常也称为序列$ \{c_t\} $与序列$ \{X_t\} $的一个\textbf{卷积}。如果我们定义函数$ C $为，
\[ C(e^{-2\pi i f})=\sum_{j=-\infty}^\infty c_je^{-2\pi if j} \]

则一些运算表明序列$\{X_t\} $与序列$ \{Y_t\} $的谱密度$ S_X(f) $和$ S_Y(f) $存在如下关系，
\begin{equation}\label{3-2}
 S_Y(f)=|C(e^{-2\pi i f})|^2S_X(f)
\end{equation}
 
\eqref{3-2}式对于AR，MA，ARMA过程的谱密度计算十分有用。
\subsection{AR，MA，ARMA过程的谱密度}
\textbf{白噪声}\hspace{1em} 根据\eqref{3-3}式，白噪声谱密度为
\[ S(f)=\sigma_e^2 \]

\textbf{MA(1)的谱密度} \hspace{1em}MA(1)就是$ c_0=1,c_1=-\theta $时（注意，此时$ j=0,j=1$，且其他$ j $处，$ c_j=0 $）白噪声的一个滤波，根据$ C $函数之定义，有，
\begin{align*}
|C(e^{-2\pi i f})|^2&=(1-\theta e^{2\pi i f})(1-\theta e^{-2\pi i f})\\
&=1+\theta^2-2\theta\cos(2\pi f)
\end{align*}

\eqref{3-2}式表明MA(1)的谱密度为，
\[ s(f)=[1+\theta^2-2\theta\cos (2\pi f)]\sigma_e^2 \]

\section{谱估计}
\textbf{非参数谱密度估计} \hspace{1em}样本谱密度不是理论谱密度的一致估计，需要其他的估计方法。回忆非参数估计的思路，可以得到一个平滑谱密度估计，
\[ \bar S(f)=\frac{1}{2m+1}\sum_{j=-m}^{m}\hat S\left(f+\frac{j}{n}\right) \]

它是对邻近的$ 2m+1 $个样本谱值取平均得到的。实际上不一定是等权重平均，还可以如下约束形式进行加权平均，
\begin{align*}
W_m(k)\ge 0\\W_m(k)=W_m(-k)\\
\sum_{K=-m}^{m}W_m(k)=1
\end{align*}
以得到，
\[ \bar S(f)=\sum_{K=-m}^{m}W_m(k)\hat S\left(f+\frac{j}{n}\right) \]

$ W_m(k) $称为\textbf{谱窗}。类似地，关键在于$ m $的选择，一般建议$ m=\sqrt{n} $，也会试用$ 2\sqrt{n},\sqrt{n},\frac{1}{2}\sqrt{n} $。在R语言中可以使用\lstinline|spectrum|
函数进行估计，里面可以选择$ m $、加权的核$ W_m $以及置信区间$ ci $。

\textbf{参数谱密度估计} \hspace{1em}实际上，可以先拟合一个时间序列模型，然后按照上一章讲的，每一个参数化的时间序列都有自己的理论谱，这就是它们的谱密度估计。